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ABSTRACT 


An investigation of the body wave travel time 
data obtained from the Apollo program was undertaken. 
The main effort was directed at finding all velocity 
depth profiles consistent with the travel time data. 
This P and S wave data sampled the first 100 km of the 
lunar interior. Using a radially symmetric velocity 
structure as an approximation to the moon this lunar 
seismic data was inverted by the Hedgehog method. 

This method was used for the investigation because 
of its ability to deal with both nonlinear inversion 
and uncertainties in the data. 

Differing assumptions on the accuracy of the 
lunar data and on some interpretations of body wave 
arrivals were made to see what effects these had on 
the acceptable range of velocity structure the Hedgehog 
algorithm produced. These were then critically examined. 
The best P wave velocity structure obtained was one 
containing the following : 

1) A Seeotwaee velocity increase in the first 4 km 
to a velocity of 4 km/sec. 

2) A continuous increase to a velocity of 4.5 km/sec 
at 14 km depth. 

3) A rise in velocity to 6.5 km/sec at 20 km 


depth. 
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A constant velocity layer of 6.5 km/sec 
eee maang to 60 km depth. 

A rise to a high velocity lid of 8.5 km/sec 
at a depth of 60 km. 

A drop to a constant velocity half space of 


7.75 km/sec below 80 km depth. 
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CHAPTER 1 
INTRODUCTION 


The basic problem of geophysics may be stated in 
the following manner: Given j observations e on the 
earth (gross earth data) and given an earthmodel m of 
n physically interesting parameters that describe a 
set of properties of the earth, and given a set of rules 
that assign to an earthmodel a set of theoretical obser- 
vations what can we determine about our earthmodel m? 
Or rephrasing this ee: Given an incomplete set of 
observations on the earth what can we say about its 
structure? | 

This is the statement of the inverse problem of 
geophysics watch this thesis will examine, especially 
with regard to the practical aspects of the inversion 
of seismic data to obtain an earthmodel. Seismic data 
is often inaccurate, inconsistent, and occasionally 
incorrect. I will examine some direct inversion tech- 
nigues and some inversion methods that rely on a search 
through an earthmodel parameter space. 

The direct inversion techniques to be examined 
shall include: 
ci Herglotz-Wiechert type inversion which involves 

the version of an exactly known travel time 

(T-A) curve to obtain a seismic velocity profile 


of the earth. 
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2) Linearization procedure due to Backus and Gilbert 
(1967-1968 a,b) which uses a linear perturbation 
theory to obtain a best fit model and then examines 
the trade off between variance and resolution in 
the model. Jackson (1972) has examined the non- 
uniqueness of the solutions for the linearization 
technique by examining the marginally acceptable 
models —- the Edgehog eee ee 
The direct techniques to be examined involve 

parameterization of a set of earthmodels. A parameter 

vector (a point in parameter space) thus defines a par- 
ticular earthmodel. The task is then to find all accep- 
table points in best Seti space with an acceptable 

point being one which produces theoretical observations 

that agree with the actual observations. Two methods 

are in oifttion use: 

Dr) Monte Carlo search which is a random search in 
parameter space. 

2) Hedgehog which is a controlled search from a 
starting point and delineates an entire connec- 
ted region of acceptable points. 

It should be noted that in seismology the obser- 
vables are the travel time T (the time it takes seismic 
energy to propagate from a source to a receiver), the 
ray parameter p = dT/dA, where A is the earth centre 
angle (Figure 1), and the amplitude of the seismic wave 


at the receiver A % a?*r/ah-. Due to the generally high 
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frequency of seismic waves with respect to velocity 
structure (velocity as a function of depth) in the 
earth the travel time may be described in terms of 
the theory of geometrical optics with exceptions 
where diffraction, interference, or tunneling become 
important. 

Applying geometrical optics we obtain the set 
of relations in Figure 1 for a seismic wave whose 
source is on the earth's surface in a radially symme- 
tric and isotropic earth. 

The theories of evolution and genesis of the 
moon are modernized equivalents of theories that fall 
into three broad categories. That is fission from the 
Carth, Capture, or formation of the.moon-from the solar 
nebula in an orbit that is in close proximity to the 
earth. 

Thevonily restrivclion that “seismic data puts, on 


these theories 16 (hat arleast jsome por tionc Of the 


moon were molten and differentiation resulted. 
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is the ray parameter; a constant along a particular 


ray; 
is radial position; 

is the velocity structure; 

the earth centre angle between source and receiver; 
the radius of the earth; 

the turning point of the ray. 


Figure l 


Seismic Ray Theory 
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CHAPTER 2 


DIRECT INVERSION 


The traditional method of inversion for a travel 


time curve is due to Herglotz (1907) and Wiechert (1910) 


with the solution (Jeffries 1962) given by: 


a (p,) 
max aL DS 
(a) T in (> ) = cosh (——) dA 
min Pi 
0 
with Py the value of the ray parameter for the ray 


bottoming out at r ). (See Figure 1 in the 


se Va) 
introduction for notation.) - This is an exact solution. 


If you know the structure above r 


min YOu can find the 


structure at rin’ Thus one works down from the surface. 

This method has an inherent limitation-in that 
the ray parameter must be known for all values of A 
between 0 and A(p,) and this method fails in the presence 
of a low velocity zone (Figure 2). 

The cases where there are low velocity zones may 
be taken care of by an extension of this method due to 
Gerver and Markushevich (1966, 1967, 1968). It is best 


described in Chapman (1973) and is summarized below. 


We know that A(p) can be written as below 


x 
0 
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Figure 2a) 


Triplication 


A triplication occurs when there is a zone of high velocity 
gradient (dv/dr>v/r) which is characterized on the travel 
time curve by three branches. 


Often these separate branches 
are so close in time as to be indistinguishable. 


Figure 2b) 
Low Velocity Zone 


When there exists a layer in the earth where SY (iar ale this 
LVZ)aiwtwo rays wi ray 
is called a low velocity zone ( ) Ae 
i ill surface with a 
ters P and p+ dP as in the figure w 
ee ee ERS. Hence a shadow zone (D-E) in the travel 
time curve where no observations (T=) exist. 
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applying the integral operator 


eS ee 
(p= Py) 
I 


to A(p) we obtain 


n r 
(4) : 2 : Bree ss. ee OD oil 
OU ied: 2 17a a ey Bg 2 
x 2 nepal tapes 
This is an integral in the p-r plane where this integral 
is over the area Da! with the area between r, and a a 


low velocity zone. For simplicity only one low velocity 


zone was shown in Figure 3a. Let 


k 


0 S* 
(S00) 2D) e007 (pa 
k=l] 
. * * 2 ° 
with Pee] < Py < Py and k* being the total number of 
low velocity zones. Reversing the order of integration 


we obtain Figure 3b. 
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where the integration is with respect to all branches 
of the travel time curve. That is in the case where 
there are several branches then the integration is 
positive on advancing branches and negative on receding 
branches. 

Examining the case where k* is zero (there are 
no low velocity zones) then we see that we have a solu- 
tion of the form 


(py) 


f A 
io y a, In (— ee | cash (by da 4... 
min 0 Py 


Knowing an exact and complete T-A curve we can obtain 
p(A) = aT/dA. If we know the structure above any Py 
we can integrate (1) and solve for Cmin’ Since 

Nt eee, = Pr/" min we then know the structure at this 
point, and using the same procedure we can find the 
structure for Py AP dp: 


For the case where k* = 1 we have one low velo- 


city zone and we must evaluate the term 


2 2 
OM 2. -1(0 hPa 

( aa ) | r tan a, dr ° 
r Py ~ Py 


Above the low velocity zone the only term to appear is 
that of equation (10). *Thus ry (the upper bound of the 


low velocity zone) can be determined. Since this is a 
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low velocity zone we can put a maximum upper bound on 
the velocity and n. From the travel time curve we 

can measure the time jump 0, across the shadow zone 
due to the low velocity zone. This is a measure of 
the extra time through the low velocity zone. Knowing 
this and the maximum velocity a least lower bound can 
be placed on ri: Thus eco ome bilan will look as in 
Pigure 3¢2 te 


Gerver and Markushevich (1967) have shown that 


through any point in this solution space (the shaded 


area) there exists one solution which goes through that 


point. Note that not all curves lying within the 

shaded region are solutions. For any curve through a 

low velocity channel the curve to the next channel is 

uniquely déveiminga, 

There are several practical difficulties in the 
-application of this method which make it unsuitable 
for inversion of travel time data. These are: 

18) There may exist branches (triplications) on the 
travel time curve. This method requires exact 
knowledge of all ray parameters on all branches 
in a given range which may be impossible since 
the different branches of a travel time curve 
are often inseparable on the seismic record. 

oy The ray parameter p = dT/dA is required. 

Unless we are dealing with arrays where a 


direct measurement of p can be made we are 
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forced to differentiate the travel time curve. 
Numerical differentiation of experimental data 
is inaccurate, This is especially important 
where the travel time curve flattens out and 


where a small uncertainty in p would lead to 


large uncertainties ina resulting earthmodel. 


This method is based on classical ray theory 
which is inadequate in regions of high velocity 
gradient due to diffraction. This diffraction 

is frequency dependent leading to dispersion in 
the wavetrain. Thus for example P and PcP data 
cannot be simultaneously inverted by this method 
to obtain the core radius. 

Uncertainties in the data are not taken into 
account. We have assumed a complete and accurate 
knowledge of the travel time curve. Unfortunately 
seismic data is a nonrandom, inaccurate sampling 
Of a travel time curve. Thus to make use of this 
method we must fit a curve to our data hoping 
that it is a near enough approximation to our 
travel time curve to calculate an accurate earth- 


model. 
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LINEAR INVERSE THEORY 


The statement of the inverse problem in geophysics 
is the following: Given J measured gross earth data 
Beg Epi ste EY with each EY a function only of an 
earthmodel (the data has been corrected for instru- 
mentation effects) m(r) = (a, (rx), ao (L)y..- Ay (2) ) what 
do we know about each a; (r)? For example the earth's 
mass, moment of inertia, travel time between two points 
of P or S waves are all gross earth data. The a,'s are 
physically interesting functions of the earth,generally 
but’ not necessarily assumed to be functions of radius 
such as density (p(r)), or P wave velocity (a(r)). 


Let us define an earthmodel m to be an ordered 


N-Tuple of functions on 0 < r <1 with 
(12) m= (CER DC) 


and define another earthmodel im' 


L3) m' = (a,, a. 


v 
slelen pice e 
Ud a a? 


Then we can define a multiplication by a scalar and 


addition as 


(14) bm + b'm! = (ba,+b'aj, ba,tb'a},...,ba,+b'a,) : 


We can also define an inner product 


i 
(15) Gy) = | delay (ral (x) a, (rag (r) +. cay (ray (r) 
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With these definitions and the realization that we can 
restrict ourselves to deal with only a finite set of 
physically interesting functions that are piecewise 
continuous and bounded then 


i | 
(16) | | |a; (4) | [dF < © 

i=1 0 
which means that we can complete this space to form a 
Hilbert space M.- | 

It should be noted that in equation (15) N 
positive weighting functions, depending on r, May be 
inserted. For example if a model is comprised of a 
series of concentric shells we may wish to weight the 
shells with a value dependent on their width. 

It should be obvious that the ith gross earth 
datum is a rule which assigns to every fhe M, a real 
number E, (im). In general then a gross earth datum 
will be a nonlinear functional on M.- Since linear 


functionals are easier to deal with than nonlinear 
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functionals consider the example of the following linear 


functionals. The mass Ey and the moment of inertia E, 


of the earth 


1. 1d 
, - 8 2 
(17) E = At | dr p(r) : E, = al Grarnep (a) < 
0 0 
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If a measurement of a particular gross earth datum 
E. for the real earth is EY then we know that the point 
mie M, which describes the real earth lies on the hyper- 


surface defined by 


>, Nee /< 
(18) E,(m) = Ey: 


Le E, (m) is a linear function then the solution 
space will be a hyperplane. 

If we have Ea Bees Bs gross earth data then 
we know that the real earth lies on the intersection in 


M of the J hypersurfaces. 


ie! A ee 
(19) Ee E. 


Thus it can be seen that solely from the measure- 
ments EY we can deduce nothing about the Spivey To be 
able to make deductions we need to know the j functionals 
E; (m). The importance of knowing the E, (m) cannot be 
overstressed (Gerver 1973) since any errors or approxi- 
Mations inherent in them will be carried throughout the 


problem. 


The Problem of Non-Uniqueness 
Since M is an infinite dimensional space and 
equation (19) introduces only a finite number of res- 


trictions intuitively we can see that our solution set 
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will be an infinite dimensional subset of M. Consider 


the following pathological example. 


1 

(20) Bl) = n+] drf(ay-1)? + (ay-2)? + 02. (ayn) 7 
| 0 

Then the observation Eo (mf) = n imposes only one 


condition yet the solution is the point 
(2h) a, = eg iis CRC ig everest vcs My 


To remove such mathematical pathologies as this, 
we define Frechet differentiability. E is Frechet 


differentiable at m if for any 6m we have 
(22) E(m+tém) = E(m) + (M,ém) + © (6m) 


where 


and M e M 


M. = QE/da. s 
‘ sE/9 5 


M is called the Frechet kernel of E at m, and it should 
be noted. that equation .(22)tis, the first terms. of. 4 
Taylor series expansion. The condition that small 
perturbations m cause small changes in observable E 

is often violated in geophysics. The travel time curve 


with a low velocity zone is such an example since small 
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perturbations in depth to low velocity zone causes a 
small shift in location of the shadow zone which 
implies that certain points on the travel time curve 
go from a finite to an infinite value no matter how 
small the perturbation. 

Let a (2) be a starting model that we guess 
and lét ioe sq (2) be the true earthmodel in that 
it furnishes all the observed data within experimental 


accuracy. Then we seek to minimize 


a3) a) (1) 98) =| Ley geen ae 
0 
sWoiee) ant a 
where ml), gh) = (ays Aoreees a.) 
= (a,tda,, atdart...,a,tda,) 
and m1) 


= (al, Aor +0078 ) aa 


To lst order we can write from equation (22) 
(24) By Gi" 4am?) = B, ee oe 


(1) 


and m mid) 


£m ff) +é6m 


I are so close that we can 


ignore ¢« in (22) then by the variational technique 


ah yO)? 
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with a: given by 


n 
gaan 


(1) 


Neglecting ¢« in (22) implies that ém may not be 
exact, but the above process can be turned into an 


iterative process with an rasa ae defined as 


j 
(25) OTD Eg Fg Pg 


If (24) were exact then 


(26) E, (mic) = E, 


and nm) would be our exact solution but as it is all 


we can say in our iterative process is that. 


(27) |E. qineD) - 2° 41 <dEs tm) - E, | 


eae < | +> (n) | 


and ém 


Thus the above technique gives an iterative method 
for finding a solution. 

Let us restrict ourselves to a set of models w 
where in equation (24) e€ is small enough to be neglected. 
Assume that Itty the true earth lies in the set of models. 


Then the natural way to assess the adequacy of a model m 
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is in terms of the resolving power. If we assume that 
at some ro the data providesan estimate of the true 
gross earth data that is a smoothed average of the 
structure around ro then for a linear functional 


a plausible and attractive smoothing is of the form 


| 1 
(28) <i (rg) > = (A (x9) sity) = | A(x xp) «ily (x) dr 
0 


and <i (rp) > is an estimate of the true earth EA 
It can be seen that 1 (x4) = Itt, (FQ) if A(r)) = 6(r-ry). 
Thus ia (rp) is a good estimate of the true earth if 
A(r,Xo) as a function of r has a peak at ro: The 
sharper the peak - nearer to a delta function - the 
better the estimate. 

Under the criterion of linearization Backus and 
Gilbert (1968) have shown that 

N 


(29) A(x) = a 


> 
b; (ro) My 
1 5 


1 


with the b, chosen such that A is wnimodular.- 


a 
(30) | dr A(r,rq) Cn Le 
0 


We will define a measure of the width of A by a 


quantity S (x9 -A) called the spread of A by: 


uk 
2 
P) S(r9,A) = 12 | (r-r9) A 
0 


2 ar ; 
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It is obvious that an A that makes a small 
spread S(rQ,A) is a good resolution function. This 
resolution function and our estimate of it (29) solve 
the problem of non-uniqueness since all models within 
the linear regime satisfying the data give the same 
estimate <i (Xp) > £or Mie) 

Associated with each estimate <M (£9) > there is 
introduced a standard error calculated from the obser- 
vational errors. It can be shown (Backus and Gilbert, 
1968) that the smaller the spread S the larger the 
standard error at any radius ro: This is a trade off 
between precision and resolution. Thus there does not 
exist a best set of b, (ro) but instead a relationship 
between the spread and error. This is the familiar 
power-spectrum trade off where you playoff the resolution 
and the precision. Whatever compromise is used in a 
particular case depends on the desired results. The 
relationship between accuracy and resolution is a func- 
tion of depth and model as well as the error estimates 
as is shown below. The spread S can be formed from a 


quadratic as below: 
(32ktens = ) : bib.S;, 


where S. is a coefficient matrix from the inner product 


of the Frechet kernels, i.e. 
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The variance (Am) 2 of the average <m(rQ)> is given 


by 


Gate once = Peer re 


where re is the covariance matrix between EY and Bi te 


can be seen that for statistically independent data C5 
is a diagonal matrix. We construct a third quadratic 
below with w being a constant 
(3:5) BJ. = 5. .cos® + C..w sino : 
13 os) 1j 

For every 6 we can find a aa which minimizes Snk 
We note the following. When @ is small Bis resembles 
S44 and we will generate delta like functions that have 
a minimal S but will have an unacceptably large (Am) *. 


As 9 increases to 7/2 S will increase and (Am) * will 


decrease. 


The Edgehog 

The variance is a statistical way of speaking of 
the extent of our model set about our "best fit" model. 
Upon examination, the variance about our best fit model 
does not really give us the boundaries to our model 


space. As was seen previously the solution is the space 
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within N intersecting hypersurfaces and thus has 
definite boundaries to within experimental accuracy. 
Let us now then try to find the boundaries to our 
solution space. That is the marginally acceptable 
solutions (Jackson 1973). 

Let us assume that we have observed n gross 
earth data EY and that we have derived j good enough 


-relations between a model m = (m).-.m_) and the 


observations such that 


Then we can define an error e, as 


RG). 26 By - E, (i) 


1 1 


which can be expanded in a lst order Taylor series as 


below 
OE. 
safe Aa bean et be i >(1) 
(37a) e205 E. E, (m ) “L (—7yy) (sm ) 
j=l om 
| 7 EA ait) ae (1) 
where we have as in equation (23) m+ 6m =m ; 


That is we have a true earth m which we bYia pLiors 
knowledge have modeled as afl) and are incorrect by 
6m. This is equation (23) rewritten with an error 


term. Which when rewritten in matrix form is 


(37) €=y-AX 
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where ¢« is the vector of the e,'s 


is the vector of the Vlg E. - E , (it) 


ake 


A is the matrix of the AS OE LAME EyaN 


Aes 


is the vector of the om. : 


If we now work with standardized data - data 
with zéro mean and a unit variance - then as was 
previously shown our model space - the set of models 
that will produce the observed gross earth datum - 
is non-uniquely described by an acceptable model with 
a variance. Thus our solution space consists of the 


models with an r.m.s. residual r given by 


(38) gots 


In choosing an acceptable model we minimize r 
for a given resolution. We have also applied an aver- 
aging process to our model by some averaging distribu- 
tion as in equation (28) which states that as long as 
we have a 6m small enough - i.e. a model not too far 
away from a true earth - the averaging over the model 
will give us gross earth data that are within the 
acceptable bounds of the observed data. For standar- 


dized data this can be written as: 


m 
(39) s=/z j mes 
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Equations (38) and (39) then correspond to conditions 
on the variance and the resolution. 

In the Edgehog procedure we seek marginally 
acceptable solutions. That is we wish to know the 
bounds of our solution space. Thus we must solve the 


equation 


r¥= } and S +} 
or 


regal and gee Ty 2 


This linear inverse theory, as applied to 
seismology, though very attractive, has one great 
practical limitation. That is the assumption that 


small perturbations in model ém cause small changes 


in observable E|- Such is the case where there are no 


low velocity zones. Even where there are no low 
velocity zones we have the added problem of numerical 
precision in our calculations. That is small enough 
perturbations ém used in a numerical calculation may 
produce an E, + 6B; with a numerical uncertainty >oE,. 
Linear inverse theory is attractive in that: 
Mh) Incomplete data may be used. Points on a 
. travel time curve rather than the curve 
itself may be inverted. 
2) Inaccuracies in the data are taken account 
of and appear as uncertainties in the models 


produced. 
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The 


"rules" used for the calculation of the 


observables from the earth model are not 


restricted to classical ray theory but can 


be 


we 


of 


we 


any physically attractive set. Thus if 
find we are dealing with a nonlinear set 
observables such as a travel time curve 


can apoly a transformation to obtain a 


linear set of observables. For example, we 


can examine the t-p rather than the T-A space 


Or 


observations (Bessonova et al 1974). 
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CHAPTER 3 


INDIRECT INVERSION 


Basically the indirect method of inversion 
involves a search through a parameter space; each 
point in this parameter space being an earthmodel. 
These parameters Py 1Pore++rPy are related to the 
earthmodel m = (a, (E) a5 (EB) p02. pay (ED) of the previous 
chapter in the following way. The ay (2) can be 
approximated by a set of coefficients and a function. 


That is 
> + 
a, (r) = END ame ne) 


where the p,'s are not functions of r. Knowing the 
parameters and the functions we then have an earth- 
model. Thus the parameters themselves describe an 
earth, and any point in the parameter space describes 
an earthmodel. For example if we are interested in 
the inversion of travel times we may choose aN layer 
earth with 2N parameters; the parameters being the 
depth to the ith layer and velocity in the ith layer. 
Any set of physically interesting and numerically 
compact quantities may make up a parameter space and 
thus a set of earthmodels. 

The search in parameter space is to find the 


set of all acceptable points in this space. An 
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acceptable point being one which produces a set of 
theoretical observables that matches within uncer- 
tainties a set of physical observations. This set 
of acceptable points in parameter space is our 
solution. 
The actual detailed process of the finding 

of acceptable parameter space points is diagrammed 
in Figure 4. The following is condensed from 


(Keilis-Borok and Yanovskaya, 1967). 


Observational Data 


This is the input of the physical pertinent 
data to a particular problem in geophysics. We wish 
to find all earthmodels that will, to within experi- 
mental uncertainty, produce theoretical observations 


to match these observations. 


A priori vlinits 

Fortunately we do not begin our search of 
parameter space without any a-priori data. That is 
we know from previous physical observation, or 
deduction, or by the physical laws themselves, that 
we have some bounds on our parameter space wherein 
lie all acceptable models. We know in what broad 


limits we can find the seismic velocity structure of 
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the earth. Thus we need not search all of parameter 
space but only a section thereof. This imposition 

of a-priori limits is very important. If we choose 
broad limits then our calculations of earthmodels 

and observables may be too costly. If we choose a- 
priori limits that are too narrow we will have excluded 
acceptable points in parameter space from our search 


and may miss physically significant earthmodels. 


Parameterization of Earthmodels 


We can approximate a ead earth by a finite 
set of parameters. The choice of parameters depends 
upon the problem and data being considered. For 
example in the case of inverting seismic data (P or 
S wave travel time, or ray parameter, or amplitude, 
or surface wave dispersion, or free oscillations etc.) 
the parameters usually chosen are some set of descrip- 
tors of 0 Ce) (P wave velocity in the nth layer), 
Bf) (S wave velocity in the nth layer), p,, (4) 
(density in the nth layer), Q, (4) (attenuation in 
the nth layer), La (location of the nth layer). Thus 
the parameters will be physical relevant descriptors 
of some set of functions which describe the earth. 

We must avoid two extremes. A too simple parameteri- 


zation wherein we would miss physically relevant 
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models and a too complicated parameterization which 


would be too costly in terms of computer time. 


Choice of an Earthmodel 


Recall that we are looking for a region, or 
set of regions, in parameter space that contain only 
acceptable points. Pihevonly practical way, to 
accomplish this is to use an N dimensional net (with 
N the dimension of the parameter space) and instead 
Of Searching An ditinity of vpoints ina region just 
to examine the points on the nodes of the nae: 

The step size of this net in any dimension is 
obviously important. Too fine a step size and the 
computations become too costly; too coarse a step 
size and we may lose information. 

It would be desirable to examine all nodes on 
this net. Generally this is inefficient. Instead 


we may proceed in one of two ways: 


Random Search - Monte-Carlo Inversion 


This involves the picking of random points in 
parameter space and has the advantage (Press 1968) 
that very little subjective bias in the choice of 
models is introduced. A good random sampling of 


parameter space bounded by a-priori limits is 
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searched, hopefully without the imposition of any 
subjective physical bias. The disadvantages of this 
method are that itis extremely costly, and that a 
random search by its very nature may have missed 


_physically interesting models. 


Controlled Search - Hedgehog Inversion 


This involves the a-priori knowledge of at 
least one good node in parameter space (this node 
might be found by Monte-Carlo search) and then 
examining in turn all the nth order nearest neighbours 
of this point. All the neighbours that are accepted 
are then examined in the same manner until the region 
of acceptance is completely Holi deated (Figure 5). 

As an example of the two dimensional parameter 
space Hedgehog search and some of the pitfalls in 
choosing a wrong grid size examine Figure 5, The heavy 
lines form a coarse grid, while the light lines form 
a fine grid and A, B are the regions of acceptance 
that we wish to find and delineate. 

We define the nth order nearest neighbours of 
a point m to be those points m+, 6m, where 6mm; is 
the step size in the ith direction, e, being +1 or 0 
and only qnqeOmLessmmofethe nhc =odec 

For example if we found the point 1 (Figure 5) 
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EXAMPLE 
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would next be examined. The point 2 is an acceptable 
point. Thus the points 6, 7, 8 would be examined and 
no more acceptable points would be found. Thus using 
the coarse grid with an order one search the region 
of acceptance that would be found would be the points 
1 and 2. Region B and the points 9, 10 would have to 
be found independently. 

Using the coarse grid with an order 2 search 
SlLoartindear. (he wpolnr, ae the pOIntS.2 p35: 45. Siainliy wey 
9, 13 would be examined. The good points 2 and 9 
would be found. From 2 no good ee would be found 
but from 9 the points 10 and 11 would be found. The 
region of acceptance that would be found would then 
include: the pointseljp2,. 0G 6L0, Ll, 12... Werwould 
consider this to be a singly connected region which 
is false. 

Using the fine grid with either an order 1 or 
2 search would produce the acceptable points l, 2, 9, 
OVS aie el 5 | 

This illustrates some of the problems inherent 
in a Hedgehog search. To summarize they are: 

L) Only connected or nearly connected regions can 
be found from a single starting model. 
2) Too coarse a grid or too low a neighbour search 


and you may miss parts of the region. 
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3) Too find a grid or too high a nearest neighbour 


search and the cost becomes prohibitive. 


Calculation of Theoretical Observations 


This is the direct problem of geophysics. That 
is given an earthmodel, a node in parameter space, 
calculate the theoretical observables. The calcula- 
tion and the theory can be whatever is physically 


attractive. 


Comparison of Theoretical to Observational Data 


We wish to know when the theoretical observations 
match our geophysical observations. This is a mathema- 
tically inadequately developed problem (V.I. Keilis- 
Borok 1970) and is unfortunately highly subjective. 

Briefly we wish to know when a set of j obser- 
vations 0; match a set of j calculations C;. We will 
say that we have an acceptable earthmodel when some 
or all, depending upon the situation, of the following 


criteria are met: 
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(3) 


2 
(4) ioe (Opecion is a, 


where the n,'s are a set of thresholds and k is a 
subset of observations (i.e. those delineating a 
Eripiication)k 

We are not restricted to this set of acceptance 
criteria, but can rather choose any physically reason- 
able set. Thus we need not use a root mean square 
measure as in (1) but could have gone to an absolute 
value measure. 

When the acceptability of a point in parameter 
space has been determined, acceptable points are 
translated to an earth structure and the Hedgehog 
search continues. 

The advantages of using the indirect method 
are that we are not hampered by our ability or in- 
ability to solve a mathematically difficult inverse 
problem with their attendant assumptions. We are 
only limited by our knowledge of the physics of the 
direct problem, which is usually adequate, and by the 
number of points in parameter space we can afford to 


search. 
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A disadvantage is that we must be careful in 
our parameterization so as not to introduce a subjec- 
tive bias. As an example of this consider the case 
where a triplication in a travel time curve is 
observed. If we have parameterized an earth in such 
a way that no regions of high velocity gradient (dv/dr) 
can exist then we will not produce a theoretical tri- 
plication. The problem can be more subtle in the case 
of the existence or non-existence of low velocity zones. 
Furthermore we do not know if our physical sciences 
background has limited us so that we may subjectively 
miss adequate models that do not in some manner meet 
with our preconceived expectations. Besides this there 
is the problem that indirect methods of inversion are 
generally costly. That is they require large amounts 
of computer time and larger amounts of time in the 


development of the computer programs. 
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CHAPTER 4 


INVERSION OF THE LUNAR SEISMIC DATA 


The Lunar Data 


The lunar ‘seismic data has been obtained through 
the passive seismic experiment (PSE), a set of seismo- 
meters deployed by the Apollo missions. Seismic energy 
was supplied by the impact of the Saturn IVB stage and 
lunar excursion module for each of these missions on 
the moon. A full description of the PSE can be found 
inv Latham et ‘aly (Clo72)e 

Separate seismic phases have been identified on 
the seismograms (Toksoz et al 1972). If we assume that 
the moon is radially homogeneous, that is seismic 
velocity is only a function of radius, then we can. invert 
the travel time curves to obtain the velocity structure 
of the upper part of the moon. The travel time data 
has been summarized in Tables 1 and 2 and graphed in 
Figure 6 (Toksoz et al 1972). The data has been supplied 
through the courtesy of Latham (1973). 

Examination of reproductions of three of the 
lunar seismograms (Toksoz et al 1972) show that the 
P wave first arrivals can be picked on these seismic 
records with an uncertainty of about 1/2 second. The 


S wave arrivals are more difficult to pick. Since I 
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had access only to scaled down published reproductions 
of some of the seismograms I assumed that the three 
seismograms were representative. If this is true, the 
P wave first arrivals of most of the records could be 
picked to withan an uncertainty of 1/2 second. ) I 
assumed that the S arrival uncertainty was about twice 
this value. 

Using a radially symmetric moon small plausible 
lateral deviations can be taken into account by includ- 


ing them in the observation uncertainty. That is the 
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moon's ellipticity, small lateral variations in velocity 


structure, and elevation of the source and receiver 
above the lunar geoid can be included in an uncertainty 
in arrival time. 

A shallow seismic velocity profile obtained from 
the Apollo 17 active seismic experiment (Kovatch and 
Watkins 1973) showed that the first few kilometers of 
the Taurus-Littrow region of the moon have a step-wise 
velocity depth structure compatible with a series of 
lava flows. This seismic profile samples both mare 
and highland regions. Plausible variations in the 
first few kilometers of the moon composed of a surface 
layer of such a series of lava flows and/or intruded 
basalts suggest that differences of up to one second 
froma standard model can and should be expected. This 


is the variation in travel time caused by doubling or 
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Table 1 


P WAVE TRAVEL TIMES 


Delta in Travel Time 
Degrees in Seconds Weights 
72o7 seu 4.0 
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On 0 20 
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* Only 1 of these 2 arrivals were used. In class B 
models, which one is used affects the half space 


velocity acceptance region by %.5 km/sec. 
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Table 2 
P WAVE TRAVEL TIMES* 
Delta in Travel Time 
Degrees in Seconds Weights 
1 oh: bO £5 
eS 220 2.0 
AO A ee 220 
apd ag po tai) 
ele. 2 54-40 2.0 
S WAVE TRAVEL TIMES * 
ie 16 30 ee 
deeS ZAS ge oe 
J pene 29.5 eZ 
3°40:7 3365 | 
5 61): 60.5 ee 
6.12 65.0 Lez 
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4 Note the high uncertainty in these points. It 


is partially due to their being read off of travel 
time curves in Toksoz et al (1972) and not from 


tabulated results. 
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halving the thickness of a layer in which P wave velo- 
cities increase from 1 to 4 km/sec. 

As is well known the P and S wave velocity 
structures are related through Poisson's ratio.” I 
will did a COnsStanct  Poisson”s ratio “equal ‘to 625 
and feel that the petrology of the moon is not 
sufficiently well established to justify the choice 
of a more detailed Poisson ratio structure. In any 
case the seismic velocities are not sensitive functions 
GE FOlrssOn.s Catlo an the "range .2Z2 to .Z28 the experi~ 
mental ranges observed in lunar rocks. I also feel 
that knowing the uncertainty in S wave arrivals as 
poorly as I do that trying to put limits on the Poisson 
ratio structure would be presumptuous. Thus the S wave 
arrivals have not been treated APaieyernlavie kt of the P 


wave arrivals. 


Choice of Inversion Method 


The limited amount of seismic data available on 
the moon from the Apollo seismometers as well as the 
relative lack of certainty about the internal structure 
of the moon suggest that an exploration of the range of 
- Lunar velocity models implied by the Apollo seismic 
data might very well be done by means of the Hedgehog 
‘technique. The quasi-linear techniques are difficult 


to apply in the case of travel-time data since incomplete 


a | 36h woo t ao ' aRexoMk aed 
Wicoley evaw Re bas. q art sok Lew st mo 


a4 


votse: anaes augnis Sannin 0 weaaa 


aN ot LBupS ok vast 2" 


ene 


dnsdasos A 


te 
3 
gy 
2 
=e, 
a 
.S 
/ 
AE 
ie 
Ag 
i 
is 
& 
rei 
= 
: = 
—— 
= 3 
a 
= 


weg -sutonage objet aceniot nici ap 

anorsoney! avid fence Sete S78 esitiootey pba ion oie 

“LI9gxS Sst SS. 0 oF S€. opiat ont oes otters # ‘noaetod 
fest cals I  , exepe waHGL, ab bovasade eopany ‘i 

a6 alevina ts svew a cud \aleawoas ne patwont “ 

RoBelod one aS. % etimis: 20 oF babys ae es = as. y¥ . 

Sys 2 add ayay -adotdamiasicg od’ biuow emusoilite © ‘ex 

q ons Ne, pheipacineehi Dennens: need r6u4 ever efewte 

| | otavinas | 


43 


sets of travel-time data appear not to be well- 
conditioned to iterative searches for best velocity 
depth functions. Although it is possible to search 
the travel-time data by quasi-linear techniques, I 
prefer to explore the topology of the space of all 
possible solutions, not just those that can be 
linearized, for the lunar seismic inverse problem 
by means of the Hedgehog technique. 

The Hedgehog technique is a method for system- 
atically determining all possible combinations of 
parameters of a class of velocity models which satisfy 
a set of observed data within given precision limits. 
I deal here exclusively with P and S travel-times. 

I do not attempt to interpret the synthetic seismo- 
grams, their character has been discussed by Nakamura 
eral (19/0), it Gs my purpose. to. detemmine upper and 
lower bounds for the velocity structure of the moon. 

I recognize that interpretations Ofisthe, synchictic 
seismograms will refine these upper and lower bounds, 
particularly in’ transition, zones... .It.1s my purpose 
here to demonstrate that fairly stringent limits can 


be derived from travel-time data alone. 


Parameterization and Lunar Models 


My Hedgehog program demands that the class of 


all possible solutions be described in terms of a 
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limited (less than 16) number of parameters. These 
parameters can be any number which is plausible, or 
possibly a number which might by some investigators 

be considered a rather implausible descriptor of 

the lunar velocity structure. The conversions from 
parameters to actual velocity structures need not be 
linear. There are no practical restrictions on the 
kind of parameterization or on the kind of conversions 
that can be carried out between lunar parameters and 
lunar structure. 

The choice of parameters for a Hedgehog program 
is governed by the need for a compact description of 
all possible solutions of the inverse problem. 
Obviously some parameters are not as well suited as 
other parameters and it is usually not immediately 
obvious which are the best parameters for a given 
inaccurate and incomplete set of data. The choice 
should be made of those parameters which are best 
determined by the available data. 

There is no immediately obvious choice of para- 
meters for the Apollo seismic data. I experimented 
with a number of sets of parameters which appear to 
me to be plausible, and chose three sets of parameters 
which gave compact description of the resulting accep- 


table lunar models. 
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hi choose three classes of models. Each class 
of models being a separate parameter space differing 
from the others in its dimension (number of parameters) 
and in the parameterization (the types of models to 
be searched). This is due to the possible different 
interpretations of the data near 1000 km (Table 3). 
Each class of models was separately inverted by the 
Hedgehog technique. 

Examination of Figure 7 will show the charac- 
teristics of the 3 classes of models. Down to region 


IV they are identical. A description of the models 


follows: 

All Models | 

Region I - First four kilometers of the moon. No 
variation allowed in v(r) which is a 
smooth curve approximating the shallow 
seismic results of Apollo 17. 

Region II - A smoothed step with the position of the 


top of the step, velocity at) the! top of 
the step, and velocity gradient, free 
parameters for the Hedgehog algorithm. 
Region III - A constant velocity layer. This velocity 
is a free parameter. 
Region IV - A smoothed step with the position of the 
step and the velocity at the top of the 


step free parameters. 
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Table 3 


P WAVE TRAVEL TIME DATA AT LARGE DELTAS 


Delta in Travel Time 
Kilometers* in Seconds Weight 
Meteoroid 
impact of 
day 134, 967 Sa 20 
1972 1026 1365 1.0 
Apollo 16 
SIVB 1099 147.0 0.4 
From Table 1 850 123.0 2.0 
1032 L550 Zin 
Lo5Sk 220 


30.3 kilometers equal 1 degree on the lunar surface. 
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Class A Models 


Region V - A constant velocity half space; the velocity 


in the half space is a free parameter. 


Class B Models 


Region V - A constant velocity half space with a 
L.V.Z. at 100 km depth; velocity in the 
half space is a free parameter. Travel- 
time data in the region 800-1200 km are 
interpreted as multiples (i.e. PP travel 


times). 


Since there exists a low velocity zone near 100 km in 
all these models first arrivals near a delta of a 
1000 km will not be direct P wave arrivals since this 
will then be in a shadow for direct P arrivals but 


will be the multiple PP arrivals. 


Class C Models 


Region V - A step down to a lower velocity half space; 
the velocity at the top of the step and 
the ivelocity of the hall (space are tree 


parameters. 


Acceptance Crireria 


Our criterion for a lunar model are that the 


following 2 conditions hold for both the P wave and 
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S wave data separately: 


(1) “ge 
0; - Ct 
(2) no more than j of the | 7 =| SNE 
tt 
where 


O is the ith observed travel time; 

C is the ith computed travel time; 

W. is the ith assigned weight; 

N is the total number of points on all branches of 
the travel time curve; | 

fe) is the rejection limit for the r.m.s. deviation; 
j,K are rejection limits such that we allow no more 

than j of our observations to deviate more than 


a time k. 


There exists some difficulty in deciding which 
observation goes with which calculated phase. The 
algorithm for choosing phases to which data is to be 
matched was developed at the University of California, 
Los Angeles in 1967 (Ketlis-Borok 1971). in the 
absence of dT /dA measurements, it consists of deter- 
mining that theoretical phase which is closest, in an 
absolute sense, to the observed phase. All theoretical 


phases are computed for a given lunar model before the 
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comparison is made, those phases which follow or precede 
other phases by less than 0.5S were discarded. Since 
the theoretical travel-time computation uses power law 
interpolation, it should be observed that the amplitude 
calculations may have spikes at critical deltas. The 
algorithm for discarding unreasonably close arrivals 
uses amplitudes. Thee checked that these amplitudes 
were calculated in ranges of A in which our calculation 
can be expected to be reasonably syalisene: 

It should be noted that the relative weights W. 
have been semi-quantitatively assigned to the observa- 
tions. That is a qualitative rather than guantitative 
analysis of the quality of arrivals at differing deltas 
and of various phases has been attempted. I have 
assigned various weights to these values Wien the 
weights dependent upon both comments and assigned 
uncertainties in Latham et al (1972), Toksoz et al 
(1972) and Latham (1973). I. have also separated the 
data into three separate tables. 

The first table consists of first and second 
arrivals with no possible argument over the identifi- 
cation of ‘the “arrivals? ° That’is*there™is ‘complete 
agreement between the various groups analyzing the data 


in the picking of the seismograms for these arrivals. 
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The second table consists of possible multiples 
such as PP or SSS and some poor quality first and 
second arrivals. These multiples are in some cases 
hard to pick and there exists some ambiguity in the 
picking of the first and second arrivals. That is 
there is some argument whether these can be observed 
on the seismograms. 

The third table consists of data near 1000 km 
including data from Table 1 and arrivals from the 
impacts of a meteoroid and the Apollo 16 Saturn IVB. 
Positions and impact times of these latter two impacts 
were calculated from data not present in these tables 
by Latham et al (1972). Thus these latter two points 
on the travel time curve have a greater uncertainty 
and are in a weak sense dependent upon a lunar model 
to about 60 km depth. 

Due to the differing acceptability of the data 
in these tables several different cases had to be con- 


sidered in the Hedgehoging of the data. 


Parameter Step Size 


I ran several groups of Hedgehogs on the data 
with both a coarse grid - one where the various 
velocity step sizes were 1 km/sec and position step 
sizes were 10 km - and a fine grid - one where the 


various step sizes were 1/2 that of the coarse grid. 
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Hedgehog Runs and Results 


I ran several Hedgehogs on each of these models. 


The results are summarized below. 


a} 


2) 


Class A Models: 
I égtata not find any acceptable models of this 
type with any reasonable set of acceptance 
limits. That is the data appears internally 
contradictory in that the data points going up 
to make the triplication (either from Table 1 
Ore 2) (suggesa ashalt Sseece: velocity of. 975 
km/sec while the travel times near deltas of 
30° suggest a velocity of 7.75 km/sec. With 
this class of model I could not reconcile these 
seemingly contradictory observations. 
Class B Models: 
These models are characterized by a poorly 
resolvable L.V.Z. near 100 km depth. With an 
L.V.2. located at thate-depth-ishave-no—direct 
P arrivals near a delta of 30°. Although I am 
in a shadow zone, multiples (PP for example) 
may be possible. Interpreting the arrivals 
near 1000 km as PP the following results were 
obtained: 
i) Using Table Imdata; and the acceptance-cri= 
terion that o=2.8, j=0, k=1.2 andgascoazse 
grained search I obtained the results in 


Figure 8. 
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i) «Using TableW? and 2™data; ovV= 6.0, 4?= 0, 
k = 2.4, and a fine search I obtained 
Figure 9, 

iii) Using Table 1,2,3 data - no acceptable 


solutions were found. 


I see that in cases i) and ii) the acceptable 
models are essentially the same. I can see that this 
is so because there are fewer but more accurate data 
points in case i) than in case ii). This shows up in 
the greater uncertainty of the location of the layer 
IV interface in case ii) due to a-larger choice of 
acceptance limits 9 and k. There exists one possible 
difference in these solutions in that there can exist 
possible differences in velocity gradient in region II. 

In these two cases I have explored the possibi- 
lity that the uncertainties (o and k) should be reduced 
to 1/2 of the listed values in cases i) and ii). We 
then found no acceptable solutions. Increasing the 
uncertainty obviously increased the number of possible 
solutions. Increasing the number of possibly badly 
identified phases (j) from zero to two caused a spec- 
tacular increase in the number of solutions. A badly 
identified phase is an arrival which deviates from its 
matching theoretical arrival by more than an acceptable 


uncertainty k (equation 2). 
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If all the travel times near 1000 km (Figure 4) 
are graphed then it is obvious why there are no accep- 
table solutions in case iii) since the data appears 
to be inconsistent. 

| I can conclude the following. IF the Table 3 
data are in error 710 sec, the models in cases i)or ii) 
ane correct.’ iF the Table 3 data are more accurate 
than this I must conclude that my model is inade- 
quate. For the case that Class B Models and Figure 8 
are accurate descriptors of the seismic structure in 
the moon the high velocities in the half space should 
be noted since such high velocities raise certain strict 
petrological constraints (Toksoz 1972). 

3) Class C Models: 

These models are characterized by a thin high 

velocity lid overlying a low velocity Zone. 

Classically this would produce a shadow zone 

from 12° to 50°. Applying wave theory in this 

classical shadow Poder suspect the following: 

i) A head wave travelling along the top of the 

high velocity lid propagating with the 
velocity of the lid. 
ii) Possible leakage of the low frequency com- 
| ponents of this headwave into the L.V.Z. 
Thus I will have a high pass filtered head- 
wave. The amplitude of the headwave will 


be highiy variable due to the nature of 
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headwaves in that small changes in either 
gradient of velocity in the lid, or curva- 
ture of the lid anterface will! cause 
drastic changes in amplitude. 

iii) A tunneled wave having tunneled into and 
cut of the L.V.Z. The tunneling amplitude 
is a function of frequency, velocity con- 
trast between top of the lid and L.V.Z., 
Width. of the, did, and. of.deltak .For a 
thin lid Vl seismic wavelength this tunnel- 
ing amplitude is relatively large. This 
tunneled phase will be effectively low pass 


fi¥tered (Figure 11). 


Unfortunately the travel=time program that we 
have in our Hedgehog package does not take nto account 
either tunneled phases or headwaves. Yet if we make 
the interpretation that the data in Figure 10 can be 
analyzed as follows: 

i) two of the three phases marked * are tunneled 

| phases (I can ignore either of the points 
at 1038 km); 


ii) the three phases marked +t are headwaves; 


‘then this is consistent with a thin high velocity ilGles 
The position of the lid (and the velocity structure 
above the lid) is that of the interface in the Class 


B Models. The velocity at the top of the lid is 
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Velocity —~> 


| ig 
radius ——~ radius —~>_° 


A wave incident at the high velocity lid with ray parameter pg will tunnel energy through 
the lid which will be reflected at r-¢.Some of this energy will resonate in the L.V.Z. the 
rest will tunnel out. 
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v8.5 km/sec (between 8.0 and 9.5 km/sec) and a constant 
velocity halfisspace cof ww7 15 .km/sec !} (7525 sto <hi75 km/sec). 
If this last class of model represents the true 
moon then a spectral analysis of seismograms in the 
range 500 km § A S$ 1500 km may show evidence of such 


a lid. 


Conclusions 


Down to a radius of 1670 km all 3 classes of 
models coincide. The data gives moderately good 
resolution of the model down to this depth and we 
are fairly certain that the seismic velocity profile 
lies somewhere in or very near to that of Figure 2. 
The width of the acceptance region is highly dependent 
upon the acceptance criterion. We cannot tell how good 
these are without a careful and critical analysis of 
the seismograms. 

Below this depth I am faced with the following 
prospects: 

1) Class A Models: 

These appear too simple. 

2) Class B Models: 

For these to be acceptable some of the data 


points near A of 1000 km must be incorrect. 
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3) Class C Models: 
These agree well with all the seismic 
observations but the petrology of a high 
velocity lid followed by a L.V.Z. on the 
moon requires special, though not necessarily 


uncommon, petrological constraints. 


Relation to Lunar Petrology 


The relation of seismic velocity structure to 
lunar petrology has been discussed in both Toksoz et 
al (1972) and Latham et al (1972). To summarize their 
conclusions for my class two and three models I have 
made the following correlations using their results. 

Region I probably consists of a series of dry 
basaltic lava flows (not a lunar soil as was thought 
prior to the Apollo 17 landing) grading into the 
second region which has velocities corresponding to 
samples of lunar basalts that have been examined. 

The third region may be a mixture of gabbros, anortho- 
sites and pyroxenites. 

In the case of Class B Models the half space 
(region, V) ,.velocity is too, high except for mineyal 
assemblages such as pure olivine. But if I invoke 
a thick olivine layer - the half space - the density 
of olivine is so high that the mass of a moon with 


such an inner structure would be completely in 
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disagreement with known value derived from astronomy. 
There have been no lunar samples returned with a 
velocity greater than 7 km/sec. This could be due 
impart’ te the*rarity*or such outcrops’ on the’ lunar 
surface. There is a similar example of the occurrence 
of dunite oh the earth. 

In the case of Class C Models the half space 
velocity could correspond to pyroxenes and/or eclo- 
gites with possibly some small amounts of olivines. 
The high velocity lid over this layer could be a thin 
high density garnet layer which Anderson (1973) 
predicts would be stable near this depth in the lunar 


interior. 


Some Constraints on Evolution and Genesis of the Moon 


The theories of evolution and genesis of the 
moon are modernized equivalents of theories that fall 
into three broad categories. That is fission from the 
Garth, capture, or formation of the moon irom the solar 
nebula in an orbit that is in close proximity to the 
earth. 

The lunar seismic data and age analysis of lunar 
samples show that ats moon was molten and differentiated 
with the last episode of surface melting occurring about 
3.1 billion years ago and the earliest known episode of 


melting occurring near 4.1 billion years (Toksoz and 
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Solomon 1972). Further analysis of the lunar samples 
show that they are silicates similar in molecular 
structure to those found on earth but lacking in water 
and being enriched in the elements Fe, Rb, K, Na but 
depletedyatinCa,y Ali apPleeUioeTh, -Ba p-Sr i (Anderson, 19:72a) . 
It should also be noted that by tracing the orbit of 
the moon back into time the moon can be shown to have 
spiralled out from the Roche limit. 

Using this entire set of data as constraints 
the three categories of lunar genesis and evolution 
theories are retained as below. 

Capture of the moon by the earth, with the 
moon originally in the nighbourhood of Mercury 
(Cameron 1972) or coming from an undetermined region 
(Alfrén and Arrhenius 1972) has been shown compatible 
with the analysis of the lunar data. 

Similarly Anderson (1972b) has shown that 
formation of the moon near the earth in an orbit 
highly inclined to be ecliptic is possible and com- 
patible with lunar data. 

Ringwood (1970) has also shown that a kind of 
fission of the moon from the earth is also compatible 
with the lunar data. He states that "during the later 
stages of accretion of the earth a massive primary 


atmosphere developed that was hot enough to evaporate 
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selectively a substantial proportion of the silicates 
that were accreting on the earth. Subsequently the 
atmosphere was dissipated and the relatively non- 
volatile silicate components were precipitated to 
form a swarm of planetesimals or moonlets, from which 


the moon accreted." 
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The entire Hedgehog Package is a superstructure 

erected about a systematic search routine in an N- 
dimensional space. This routine was first used 
successfully in a geophysical context by V.P. Valus 
(1968) in Moscow, and was subsequently used as a basis 
for the present package which evolved at U.C.L.A. under 
the guidance of L. Knopoff and V.I. Keilis-Borok. The 
present version incorporating a number of small 
modifications is based on the version KBPRGI implemented 
at Edmonton by E. Nyland, and E. Roebroek. The 
description that follows is hased on the program KBCAM 
currently available on the IBM 370 at the University of 


Cambridge. 


Some features of the current University of Alberta 
version are not in KBCAM, and vice-versa. The user 
should beware. A quick look at the source will be 


profitable. 


The University of Alberta version has a capacity for 
interactive gqraphics using routines CUCGRF, GRID, 
hardware, and software. For more details refer to the 
CUCGRF writeup, and "Computer Graphics! for Application 
Programmers", 


1.1 THE HEDGEHOG TECHNIQUE 


This section attempts to describe briefly the 
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theoretical background to the Hedgehog technique- this is 
described in much more detail in section 5, and anyone 
proposing to use this package seriously should read that 
section carefully. If we consider an Earth model, we may 
define this in terms of a number of parameters - e.g, we 
may choose seismic velocities, densities, Q, electrical, 
and thermal conductivities, etc. The particular choice 
of parameters taken depends on the particular type of 
problems being considered. In this Hedgehog scheme we 
look for the solution of an inverse problem not merely as 
a set of velocity-depth, and density-depth functions but 
as a set of suitable Earth parameters. Those parameters 
which may be the most useful in a given inverse problem 
are not necessarily the obvious ones mentioned above - 


some combination of these may be more appropriate. 


It should be emphasized that the fundamental 
elements of a Hedgehog solution are parameters whose 
choice is up to the user. It is these parameters which 
are Sane nubeaes generate a suite of structures, and the 


method of perturbance gives rise to the Hedgehog. 


Within the Hedgehog method we allow a limited number 
(N<16) of the parameters to vary, and search for those 
combinations of sateutters which give an Earth structure 
whose properties agree with the observations within the 


required precision. We have thus a problem, how do we 
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generate an Earth structure from a given set of 
parameters. This is acheived within this version by the 


use of the symbolic language PEAR. 


As an example of the possible problems note that 
alpha can be derived from the density by an equation of 
state, and is related to beta by poissons ratio which may 


be a function of depth. 


For each of the N variable parameters we specify an upper 
and a lower bound and a step size. This defines a region 
of search in parameter space spanned by a N-dimensional 
network in which each node is associated with a different 
Earth structure. By some means, e.g. guess-work or 
Monte-Carlo methods, we choose on node as a starting 
point for the Hedgehog. This node may be "good' 

(i.e. the corresponding Earth structure fits the data) 
or, it may be ‘bad'. The method is to search all nodes 
in the vicinity of the initial node and to determine the 
region in the N-dimensional space which contains good 
nodes. We do this by changing each of m<N parameters one 
step in the positive direction and one step in the 
negative direction across the network until all the 
neighbours of this first base point have been 
investigated. Since we move only one step on any 
parameter the eae the Hedgehog, m, is the numbers of 


parameters varied simultaneously. Now, hopefully, some 
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of the neighbours are ‘good nodes. These we record by 
their node indices and when all the neighbours are tested 
we use the next good node as a base (nb an order of 1 may 
miss a point). The process is illustrated by figure 5 
Why the name Hedgehog ? - well the pattern of movements 
in parameter space is sufficiently spiny to suggest to 
some the shape of a Hedgehog. It is possible that the 
zone of good points is not simply connected or consists 
of distinct regions - the latter can be handled by using 
the combination of Monte-Carlo, and Hedgehog method 
available in this form of the Package. In the former 
case apparent fragmentation of the volume of ‘good! nodes 
may occur if the steplength is not small enough. 


1.2 USES OF THE PACKAGE 


BY euitanne choice of Program Procedure and 

calculation limits the package may be used for the 
following purposes: 
a) Inversion package for Love and Rayleigh wave 
dispersion 

and Q and SH travel times - these may be inverted 


separately or jointly 


B) Dispersion package - calculation of Love and 
Rayleigh wave 


dispersion and Q 


a ar) vt haboas 
astaamavod a wsedindd, Bas ‘Shee 


ht fats setae af nosepnen, rs ot ae 
7 hala =) ae barpandon ehame. ton at atnbea — — 


ag rvaety 


aes 409 fera sd vee sendory sat | lea oh nod enue a 


ib 


Lo 


boi ea vedaghen up {swllneo-aseptt Re woke 


5 ato cll ix.  29aaKt eae 2) wot ony an 


wero Be na ng 


_Rleipeitiew ae | eta . ta . 


pabinsit oi ABD aaa ed seman 


vt 


thoon a ae anuloy set: te: ott we ammyiat' 
Abie ae tiginss Oe tite arpestgeae whe ab 
neadoie how: toe 


= ay 


Sae ce a ieeaiele spbeenee 56 ‘epteds 


a 
Witoons 


= 


evew tpt oteat Sue sient 383 Jepioen a 7 7 


i! 


i bs iets 
van i er. Fe fa ee 


1 


74 


Cy Sh and P travel times - calculation of theoretical 
travel 


times for various phases as well as ray parameters 


The following notes should enable one to use the 
package but serious users should make themselves familiar 


with the source program. 


2 ORGANIZATION OF THE PACKAGE 


The basic structure of the Package breaks down into 
6 component parts 
a) Input and Control section 
b) Choice of a parametrised Model 
c) Calculation of an Earth structure form the parameters 
d) Surface wave routines 
e) Travel time routines 


f) Acceptance routine 


The package is well documented through comment cards 
throughout particularly in the introduction, which is 
reproduced for reference in Appendix 1. It is 


recommended that this should be read before any attempt 
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is made to use the package. The following descriptions 


Should also be supplemented by reading the source. 


Ze} 


INPUT AND CONTROL SECTION 


This consists of the routines 
BLOCKDATA 
MAIN 


INPUT 


BLOCK DATA: This initialises sizes of various arrays 


used in 


(Grafic) 


the 


ATA 


CSET FIONN SO ATO OOO. 


OoOMYN WLS Wh 


the package and set up default options on 


handling of input and output 


The function of the main program is to act as 


branching control to the various facilities in 


package. The possibilities are listed below. 


SELECT PATH ACCORDING TO INDEX 


BRANCH TABLE 


TERMINATE PROCESSING OF ONE BATCH AND GET NEW 


HEDGEHOG PROCEDURE 

COMPUTE CROSS SECTION. CHECK AGAINST LIMITS 
COMPUTE LOVE DISPERSION CURVE 

COMPUTE RAYLEIGH DISPERSION CURVE 

COUNT ITERATIONS. TERMINATE IF TOO MANY. 
ACCEPTS POLNT 

REJECT POINT 

SET UP PROGRAM PARAMETERS 
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10 CONVERT RADIAL STRUCTURE TO LOVE LAYER STRUCTURE 
11 CONVERT RADIAL STRUCTURE TO LOVE AND RAYLEIGH 
LAYER STRUCTURE 

2 CONVERT RADIAL STRUCTURE RAYLEIGH LAYER 
TRUCTURE 

13 COMPUTE TRAVEL TIME FOR Sh 

14 RESTART PROGRAM 


EXCUTE CONSOLE DISPLAYS 


AANHAHADaANARAAANM 
man 
Oo 


16 * MONTE-CARLO PROCEDURE 
Ls * COMPUTE LOVE WAVE DERIVATIVES 
18 * COMPUTE RAYLEIGH WAVE DERIVATIVES 
+9 * COMPUTE Q FOR LOVE WAVES 
20 * COMPUTE Q FOR RAYLEIGH WAVES 
Zt * COMPUTE LOVE WAVE DISPERSION CURVE BY 
DERIVATIVES 
See * COMPUTE RAYLEIGH DISPERSION CURVE BY 
DERIVATIVES 
€ 
e.1 29 + CALCULATE GRAVITY FIT 
30 + COMPUTE TRAVEL TIME FOR P 
Cc 


The control of the branching is acheived by the contents 
of the array CNTL which is set by the program procedure. 
In section 3 we describe how to set up such a Program 


Procedure for the package. 


INPUT: As it's name suggests this routine is the 
general 
input facility to the package. Control within 
the input 
routine is acheived by a transfer table listed 
below 
Cc | 
C BRANCH TO READ GROUP 
* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 
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C TRANSFER TABLE 

Cc 

C L = 1 PARAMETERS 

Cc 2 CONLIM. LIMIT PARAMETERS FOR VARIOUS 
CALCULATIONS 

Cc 3 CROSS SECTION DATA 

¢ 4 CLEAR PROGRAM 

c 5 LOVE DATA 

¢ 6 RAYLEIGH DATA 

Cc 7 BODY WAVE STATION DATA 

C S ARRIVAL TIMES 

Cc 8 PROGRAM PROCEDURE 

¢ 9 GRADIENT LIMITS 

Cc 10 END OF DATA 

€ 11 PNTSBL RECOVER GOOD POINTS ALREADY FOUND 
C 12 ‘STOEL TERMINATE EXECUTION OF THE PROGRAM 
c 13 SSYSIN ALTER INPUT DATA SOURCES 

C 14 SYSOUT ALTER OUTPUT DATA GROUPS 

Cc 15 IOALG CHANGE READ-IN ALGORITHM FOR DATA 
Cc 16 OUTFLG CONTROL OUTPUT OF DATA 

Cc 17 * SURFACE WAVE Q DATA 

Cc 18 WRITE COMMENT CARD 

Cc 

Cc 19 LOAD P TRAVEL TIME DATA 

Cc 20 + GRAVITY DATA LOAD 

c 21 ICNLIM PARAMETERS READ IN 

Cc 


The control of the branching here is acheived through the 
first four letters of a control word which are compared 
with values held in the array KEYTBL. These are in fact 
the first four letters of the descriptions above except 
for 17: QDAT/<18: COMM; “and 19:)P BOs. By this means the 
required data may be read into the appropriate groups of 


variables. 
see CHOICE OF A PARAMETRISED MODEL 


* A feature of KBCAM not implemented at University of 
Alberta 
+ A feature under development 
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This consists of two routines 


*,BCSET 
HGEHOG 
* NCSET: This is a routine to choose points and nodes 
in parameter space by the use of the Monte-Carlo 
method. 


There are five possible schemes determined by 


the 
setting of the flag MCALG. 
MCALG = 1 Monte - Carlo a good point 
= 2 Monte - Cario for a good node 
= ae Monte - Carlo for a good point then move to 
nearest node and use the Hedgehog technique 
= 4 Monte - Carlo for a good point and then to 
Hedgehog 
only if the nearest node is good 
= 5 Monte - Carlo for a good node then move to 
Hedgehog 
except in cases 1,2 it will be noted that the 
routine 
is intended to choose a starting point for the 
Hedgehog technique to enable searching of 
multiple 


regions to be acheived a list is built into 


HEGEHOG 
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( entry NODCHK ) to see if a Monte - Carlo 
node has 


been looked at before. 


HGEHOG: This is the basic Hedgehog routine which 
generates a 
sequence of nodes in the parameter space to be 
used to 
produce Earth structures. This routine has 
three entry 


points HGEHOG, PNTPCK, and RESTRT. 


PNTPCK: This stores acceptable nodes in an array. 
RESTRT: This can be used to enter a set of acceptable 
node for restarting an unfinished run. 


2-3 CALCULATION AN EARTH STRUCTURE 


This also consists of two routines 
CRSCHK 


LVCVT 


CRSCHK: This routine uses the input from the symbolic 
language PEAR (see section 4), and constructs a 
cross 
section i.e. an Earth structure corresponding to 
* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 
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given set of parameters, and checks that this 

structure lies within the limits (upper and 
lower 

bounds, and gradient limits) specified in the 


input 


LVCVT: This is a conversion routine using Earth 
flattening 
techniques to convert radial Earth models to Love 
and 
Rayleigh flat layer models. For Love waves the 
exact 
method of Biswas and Knopoff is used for Rayleigh 
waves 
the empirical method of Bolt and Dorman is used 
2.4 SURFACE WAVE ROUTINES 
There are five routines in this section 
TSW 
SWEVAL 
SWRAYL 
SWLOVE 


* SWDER 


TSW: This acts as a control monitor on the handling 
of the 


othersurface wave routines which are accessed from 
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MAIN 


through TSW. 


SWEVAL: This tests a theoretical dispersion curve 
against a 


data set, and produces a flag RTCNTL. 


RTCNTL = 0 if structure is good 
= 1 if structure is bad 
SWRAYL: Fast Rayleigh wave dispersion progran 


- Knopoff's method 


SWLOVE: Fast Love wave dispersion program 


- Thomson-Haskell method 


* SWDER? This is a multipurpose routine which 
calculates the 
structure derivatives of the ers curves 
far <a 
Particular model - these may then be used 
either to 
calculate dispersion curves for close 
Structtires opeto 


calculate surface wave 1/Q by Anderson's method 


Fach of the routines SWRAYL, SWLOVE, SWDER calls SWEVAL 


* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 
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so the the RTCNTL flag is set after entering any of these 


routines. 


2.5 TRAVEL TIME ROUTINES 


This section has ten routines 


TST fre 
DELFS DELFP 
DELRS DELRP 
DCORRS DCORRP 
TSEVAL TPEVAL 
PoTrfrrrs These routines calculate the travel times, 
ray 
parameters, and amplitudes for each phase 
arriving 
a given distance. Mode conversion is neglected 
so 
that these are strictly for SH/P waves. These 
routines 
call TSEVAL/TPEVAL, and so set a RTCNTL flag. 
DELFS/DELFP: These calculate delta given the ray 
parameter. 
DELRS/DELRP: 


DCORRS/DCORRP: These provide a depth correction to 


delta for 
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a given ray parameter. 


TSEVAL/TPEVAL: These complicated routines compare 
theoretical, 
and observed phases at each station, and 
determine whether the fit to the data is 
good, 


producing a flag RICNTL. 


WW 
i) 


RTCNTL if structure is good 


= 1 if structure is bad 
Note the use of ICNLIM (see the comments in BLKDTA). 


2-6 ACCEPTANCE ROUTINE 
ACCEPT: This indicates the production of a suitable 


parameterised earth structure. 


2.7 + GRVTST 
This calculates the gravity profile from a density 


model. 
a THE PROGRAM PROCEDURE 
3.1 DESCRIPTION OF THE PROGRAM PROCEDURE 


The program procedure which controls the sequence of 
operations in the running of the package consists of a 
* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 


; f hy rie oh : ; i ' Aa i ‘ Ma 
rey bowels hedaman, sete et aa 
hades om - Ot ae et ' ve Ae | 
mt oy : Cp han r ir f oi iw | ia 


baa: ma 8: isle +8. esendy borreae. she i ar . 
at mee Pre oF +P eae saan oatiesdses eo 


5 ~~~ sh 


apis, j : i 


= PEND. psi? & ) pd Pouberg :™) D#y 
‘wenn: on. 


Dot ar sah ran THe. 4} a 


i 
5 


bad rE detitalisse “es 5 / bein sal 7 Tee: 


+ (APIA ES vad atheros sug aie Ladle a a oan 


a ee) 


a at mom 
sidedhes 8 90, dottombeae odd, soya ink abet 


SuNsSeItA Atieq: feptiea A 


v4 


Ytiaish & wow etbiorg Nose one soselvoteg ser 


a ae ee 
ssouncon9 nenniene heads 


‘wane r30ma mason swe =o: wor raves 


to sogéenpee oo sor dia dotaw wautievong etn 


| 
_# to canine wpading odd ae a aa 
ao USieeeviad 2 tegne ep ign tom pes te aN ko oe 


84 


sequence of integers, each having a definite location in 
the sequence e.g. 

Sny Ser 6cn 7001 467400 42950 1072 22. 
Locn. 1 eta 3 FAL 6 7 BLEGe «28 


The instruction integers j have the following 


Significance 
j< 100 - Branch to label j*1000 within Main 
i> 700 - Jump to location (j-100) in the progran 


procedure and perform the instruction located 
there 
Following a branch to a routine which calls upon a 
evaluation, and therefore sets an RTCNTL flag. There are 
two conditional instructions depending on whether the 
test was good or bad. For example consider 
4 108 128 
4 - Love wave dispersion calculation Zs testing 
108 - RTCNTL=0, structure good GO TO INSTRUCTION 8 
128 - RTCNTL=1, structure bad GO TO INSTRUCTION 28 


The branching table for the instructions j is 


SELECT PATH ACCORDING TO INDEX 
BRANCH TABLE 


1 TERMINATE PROCESSING OF ONE BATCH AND GET NEW 
ATA : 
HEDGEHOG PROCEDURE 
COMPUTE CROSS SECTION. CHECK AGAINST LIMITS 
COMPUTE LOVE DISPERSION CURVE 
COMPUTE RAYLEIGH DISPERSION CURVE 
COUNT ITERATIONS. TERMINATE IF TOO MANY. 
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7 ACCEPT POINT 

8 REJECT POINT 

9 SET UP PROGRAM PARAMETERS 
0 

1 


Cc 

€ 

Cc 

C54 CONVERT RADIAL STRUCTURE TO LOVE LAYER STRUCTURE 
curt CONVERT RADIAL STRUCTURE TO LOVE AND RAYLEIGH 
Cc LAYER STRUCTURE 

Cy V2 CONVERTS RADIAL STRUCTURE RAYLEIGH LAYER 
STRUCTURE 

eB 13 COMPUTES TRAVEL TIME 

c 14 RESTARTS PROGRAM 

G7 15 EXECUTE CONSOLE DISPLAYS 

Cc 

ce 16 * MONTE-CARLO PROCEDURE 

ot uy * COMPUTE LOVE WAVE DERIVATIVES 

Strive * COMPUTE RAYLEIGH WAVE DERIVATIVES 

Cc; 19 * COMPUTE Q FOR LOVE WAVES 

cP 20 * COMPUTE Q FOR RAYLEIGH WAVES 

Gee t * COMPUTE LOVE WAVE DISPERSION CURVE BY 
DERIVATIVES 

S22 * COMPUTE RAYLEIGH DISPERSION CURVE BY 
DERIVATIVES 

C 

¢ 

Cry 29 + CALCULATE GRAVITY FIT 

Gs 30 + COMPUTE TRAVEL TIME FOR P 

c 


3.2 EXAMPLE OF A PROCEDURE 

As an example of a Program Procedure we give below and 
annotated version of a procedure which uses many of the 
features of the package (We will use the shorthand "In" 


for the instruction n) 


INSTRUCTION 
1 9 set up program parameters 
zZ 3 COMPUTE A CROSS-SECTION FROM PARAMETERS 
6: 6 COUNT NUMBER OF ITERATIONS, TERMINATE IF 
TOO MANY 


* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 
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CONVERT A RADIAL TO A LOVE WAVE STRUCTURE 


COMPUTE LOVE WAVE DISPERSION CURVE AND 


AGAINST DATA 
RRCNTEL =O, GO°TO £8 
RTCNTL = 1, GO TO 126 ‘reject! 


CONVERT A RADIAL TO RAYLEIGH WAVE 


COMPUTE AND CHECK RAYLEIGH WAVE DISPERSION 


i 
© 


RTCNTL GO TO 112 

RTCNTL = 1 GO TO 126 ‘reject! 

TRAVEL TIME ROUTINE 

RTCNTL = 0 GO TO 130 'provisional accept! 


RYCGCNTE = 1- GO TO £26 "reject ' 


COMPUTE CROSS-SECTION FROM PARAMETERS 
COUNT ITERATIONS TERMINATE IF TOO MANY 


CONVERT RADIAL TO LOVE STRUCTURE 


COMPUTE LOVE DISPERSION BY DERIVATIVES 


AND CHECK AGAINST DATA 


OF RE6OLTOVEZT 


i 


RTCNTL 


i 


RTCNTL 7. 5GO TOULZ6 reject * 
CONVERT RADIAL TO RAYLEIGH STRUCTURE 


COMPUTE RAYLEIGH WAVE DISPERSION BY 


AND CHECK AGAINST DATA 


“ny onodane ee in Z 


oh) ite toate | 


: p rinbaay 


nat ant ee me on an Soa E 

vie noesaTAR, ot Peace, s mate 

ier aIatd Use BOE bie! aK apveasa. 
a De ira os Reon 
gets + aer om oe Wp ts aT 
shape. ante seer 


ft TALS 
ae 7 


ome fennbe iworge Gey om i gi ports ai t 
tne pee aha, ot ao Dee = See 
; we a aay 


Age ti 


a | : ty 
Va ye agek ; F 
2 ‘ee Leen f 
a) C “ 
eu rs Ea in al : 
} a q : 


2 ie 


tarde 18, eae vos eviews o 

: Tereat Renraon AORN awa My a 

+ a alli HSH Einog a fesintia | 
"ta ant aN ee ‘b= e008 


o L 


ands: oh te merg zat or dae, Pos 
fe-la biden | aah ‘avae eae a, 


ve q ATAG  ReyTWOA) | k : * eee nt: ; ) 
a . bt CP - hae xt i : 
! ; * : Toa oe 


23 


24 


38 


39 


12 


RTCNTL 


RTCNTL 


reject 
REJECT RO 
OBTAIN NE 
GO TO I2 

derivat 


CALCULATE 


RTCNTL = 


RTCNTL 


CALCULATE 


RTCNTL 


RTCNTL 


accept 


ACCEPT RO 


OBTAIN NE 
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0 GO TO 112 'travel times! 
4+\GO.TQ L126 ‘reject! 

sequence 

UTINE 

W PARAMETERS FROM FROM HEDGEHOG 
're-iterate'! 

ive calculation 


LOVE WAVE DERIVATIVES FOR A GOOD 


Oi CO PO) 133 
1 GO TO I26 ‘reject! 


RAYLEIGH DERIVATIVES FOR A GOOD 


0 GO TO 137 taccept as a good 


141,GO PQyuZ26q4regjgeet! 


routine 
UTINE 
W PARAMETERS FROM HEDGEHOG 


're-iterate' 


The purpose of this procedure is that having once 


found a good structure, computing time should be saved by 


using derivatives calculated for this structure to look 


at nearby structures. A 


S soon as a ‘bad! structure is 
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found we return to full calculation as written, if during 
calculation of the derivatives a bad structure flag is 
obtained the previous ‘good! point is rejected on the 


grounds that it is only acceptable 


4 SETTING UP A STRUCTURE USING PEAR 


The purpose of the coding language PEAR 
(Parameterised EARth models) is to enable a user to 
specify a model in terms of parameters he wishes and then 
to be able to deduce Earth structures of a form 
acceptable for the evaluation routines in the rest of the 


package. 


The way this is acheived is to introduce a set of 
256 notional registers into the package, and to use these 
to compute from the parameter array (XCHH), generated by 
the Monte - Carlo or Hedgehog procedures, values to place 
in the Earth structure array (EARTH). The operations on 


these registers are carried out using PEAR code. 


PEAR is probably more than adequate to introduce any 


plausible constraints into the Inversion procedure. 
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Particularly one may specify a relation between alpha and 


beta which depends on the layer number, one can introduce 


an equation of state such as the Murnaghan equation or, 


One may study quite complicated "parameters" such as the 


surface wave impedance. 


4.1 THE LANGUAGE PEAR 


There are 29 


4 letter mnemonic 


are: 

CLAE, 
ADDE, 
SUBE, 
MLTE, 
DIVE, 
EXPE, 
STRE, 
IPTE, 


END, 


CLAH, 
ADDH, 
SUBH, 
MLTH, 
DIVH, 
EXPH, 
SLRE 


PETE 


CAG; 
ADDC, 
SUBC, 
NLTC, 
DIVC, 


EXPC, 


instructions in this code, each being a 


for a certain set of operations. These 


CLAR 
ADDR 
SUBR 
MLTR 
DIVER 


EXPR 


Each instruction card has the following format 


AY 


INSTR 


—. F0 


| J 


BW SE Et oh CF VOSS  f sPtes 5 


N1 | N2 | T 4 T2 


The various fields are handled as follows 


INSTR: 


The coded instruction INSTR is always present - 


the number of the others depends on the operation 
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J: Associated with this system we use a set of 256 
registers and the index J always refers to the current 


register 


Each instruction operates in terms of a Variable Var 


defined in terms of N11, N2, T1, T2 as follows 


Class E commands (XXXE): 
Here N1, N2 are the indices of the 
location in the array EARTH 


VAR = EARTH(N1,N2) 


Class H commands (XXXH): 
Here N1 is the index of the location 
in the array XCHH 


VAR = XCHH(N1) 


Class C commands (XXXC): 


VAR = T1 


Class R commands (XXXR): 


Here VAR is the contents of register N1 
For each of the classes there are 5 types of command 


CLAN?’ *Cillea PY se 91 ster *d and add Var 
ADDX: Add VAR to the contents of register J (CJ) 
SUBX: Subtract VAR from CJ 


MLTX: Multiply CJ by VAR 
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DIVX: Divide CJ by VAR 


EXPX: Exponentiate CJ by VAR 


There are in addition 4 assignment instructions 


STRE: Store CJ in EARTH(N1,N2) 

STRT: Store CJ in EARTH(N1,N2) but raise an 
emror flag Tf T1 > 3d, or 82°< CJ 

IPTEs Put T1 into EARTH(N1,N2) 


IPTRS Put T1 into register J 


The processing of the PEAR instructions is terminated by 


the instruction END. 


4.2 AN EXAMPLE OF THE USE OF PEAR 


In the subsection CROSS-SECTION DATA of the INPUT 
data which contains the PEAR instructions, two: further 


pieces of information have to be supplied: 


DTYPE = 1 Radial structure 

= 2 Love wave structure 

= 3 Rayleigh wave structure 
DNUMBR - the number of layers. 


The Love and Rayleigh wave structure are flat layer 


structures while the Radial structure consists of 


Spherical shelis. 


For reference we give the locations in EARTH (N1,N2) 
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DTYPE 1 2 3 


N2 
Thickness 1 6 aie 
(radius) 
S velocity Z 7 12 
P velocity 3 8 13 
Density 4 9 14 
Q 5 10 15 


The location EARTH(J,5) contains QBETA for body waves; 
while EARTH(J,10) and EARTH(J,15) refer to QBETA AND 


QALPHA in the Jth layer for surface wave Q. 


EXAMPLE: 


CROSS SECTION DATA 


1 16 

LPLE 1 1 6371.000 
LPT 2 1 6356.000 
IPTE 2 4 6356-080 
LRU 4 1 6336.000 
LBTE 5 1 6336.000 
IPTE 12 1 5700.000 
IPTE 13 1 5610.000 
IPTE 7 14 1 5475.000 
IPTE 15 hk 5375.000 
TEES 16 1 5150.000 
CLAE 1 5 1 

SUBH 1 1 

STRT 1 6 $ 5721.90 

637T a0 

STRE 1 1; 1 

SUBH 1 y: 

SERL 1 ) 8 1 SEZ AQ 

6371.0 

STRE We 9 1 

SUBH 1 3 

STRT 1 10 1 SPZ1.10 


6377.0 
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LPTE 
a by 
IPTE 
LPR 
TeTR 
CLAE 
MLTR 
ADDR 
STRE 
STRE 
CLAE 
MLTR 
ADDR 
STRE 
STRE 
CLAE 
MLTR 
ADDR 
STRE 
STRE 
CLAE 
MLTR 
ADDR 
STRE 
STRE 
IPTE 
SPTE 
TPTsE 
Bw i Bo 
DEPE 
IPTE 
Bd a 
IPTE 
LETS 
TREE 
TPTE 
Lees 
IPTE 
Lees 
EPRDS 
IPTE 
END 
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100.000 
100.000 
100.000 
100.000 
100.000 
80.000 
80.000 
80.000 
80.000 
400.000 
800.000 
800.000 
800.000 


2000.000 
2000.000 
2500.000 
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a) The use of H commands to generate interpolated layers 


in 


structures as well as tests on structure through STRT 
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b) the use of equations of state to determine P wave 
velocities from S wave velocities and densities from P 
wave velocities 

c) the straight forward setting up of unperturbed Earth 


structure factors. 


s) THE THEORY OF THE HEDGHOG TECHNIQUE 


This section describes more fully the concept of a 
cross-section and the use of parameters to describe Farth 
models and also the Hedgehog technique itself. The 


description is based on Keilis-Borok(1972) . 


It is practically impossible to represent any of the 
functions alpha, beta, density, Q analytically in the 
whole depth interval, and in any case the functions may 
be discontinuous. We therefore have to divide the 
structure into some depth intervals (layers), and 
separately approximate each physical function in each 


layer, by the same type of function if possible. The 
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parameters of the structure will be the parameters of the 
approximating functions as well as the depth of the layer 
boundaries. 

Parametrisation must be correct, i.e. correlated with the 
physical task and with the pecularities of the 
observational data being used. If the object is to find 
out whether some elements of the structure exist ({low- 
velocity layer, or a discontinuity of velocity within 
some depth interval, etc,), then the assumed system of 
parameters has to allow structures both with these 
pecularities and without them. The probability of 
acheiving the desired pecularity by random methods must 


be about a half. 


Clearly for correct and optimal parametrisation it 
is important to investigate as a start which features of 
the observational data are connected with the different 


properties of the structure. 


Any measure of the discrepancy between observations 
and the computed properties of seismic waves depend on 
the parameters of the structure, because the computed 
properties depend on them. Our solution is a set of 
combinations of the unknown parameters (i.e. a volume 
in the space of the unknown parameters) for which the 


discrepancy is sufficiently small. 
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The simplest and the only absolutely reliable method 
is to divide the investigated multidimensional region by 
a net, to calculate the function at each node of the net 
and to choose that point in which the function is 
sufficiently small. The simplest method of construction 
of the net is to place the nodes along coordinate axes at 
equal intervais. The step of the net must have a 
magnitude of the order of the error that is allowed for 
in the determination of the boundaries of our minimun 


region. 


If the number of parameters is large the number of 
nodes becomes so great that the described method is 
practically impossible. Thus some form of guided search 
method becomes attractive and the most generally useful 


has been found to be the "Hedgehog! method. 


A cross section is defined as the set of parameter values 


describing in some way an Earth structure. 


The Hedgehog technique procedes as follows: a single 
point of the minimum region Y { Y(1),-..--,Y¥(N) } - the 
Y(i)'S being the parameters of the cross section As found 
by some means e.g. the Monte-Carlo technique; and then 
the neighbouring points Y(j) + sum over (j) alpha (j) dY (3) 
are tried where alpha (4) = 0 or + 1 and different 


combinations of the alpha(j) are tried in turn, the dy 
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above are the step lengths on the net assumed to be 
embedded in the parameter space. Points falling within 
this minimum region are then selected and the same 
procedure applied to these points until the whole region 
is covered. When this has been done we return to the 
Monte-Carlo technique (omitting of course the region just 
found from further investigation) and look for another 


minimum region, and so on. 


The most reliable approach is to test all 
combinations of the alpha(j) i.e. to test all the 
neighbouring points which can be obtained from an initial 
point by not more than one step along each coordinate 
axis. This is however is computationally difficult for 
farcesgnuea.cach initial point having (3(N)-1) such 
neighbours. For this reason the method includes the 
possibility of making the limitation that only those 
neighbouring points for which not more than m (<N) of the 
alpha(j) are non zero are tested. lanes only those points 
which may be reached by taking at step along not more 


than m coordinate axes from the initial point are tested. 


FIGURE 5 illustrates the Hedgehog method for a two- 
dimensional space. Heavy lines represent a grid for 
which the horizontal increment is twice that for the thin 
lines. The area of the minima consists of two closed 


parts A and B. Suppose that we found point 1 by the 
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Monte-Carlo technique, assume m = 1 3; then we should try 
the neighbouring points 2, 3, 4, 5.. Point 4 belongs to 
our minimum region, so we then try points 6, 7, 8 around 


2c. 


To find the rest of our minima we have to return to 
a random search until we find the points 9, or 10 and 11 
or 12. However we would think of points 1, 2 and 9, 10 
as belonging to isolated regions; until we check by 
taking finer steps along the dotted path or introduce the 
finer grid. For m = N = 2 we would however try at once 
the points 2,3,4,5,6,7,9,13 around point 1, and then 
points 10,11 would be found from 9; in this case we would 
think of A,B as a simple area until a finer grid was 


introduced. 


This example shows that if the increments dy(j) are 
too big or m is too small, we can, easily miss some 
points of the area of minima or introduce a fictitious 
division of this area into isolated parts. However 
decreasing dy(j) and increasing m leads to a large 
increase in the number of computations. Decreasing dvy(j) 
is particularly inconvenient: it gives a lot of internal 
points in the area of minima which are of no interest. 
It is better to increase m or preferably to rotate the 


cordinate axes. 
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6 THE CONSTRUCTION OF A DATA SET 


We describe here a typical example which gives an 
indication of what is required. An examination of the 
test routines in KBCAM.DATA.TEST will give the details of 


Suitable data sets for most purposes. 


COMMENT: Title of data set etc. 
CLEAR: Initialising array values 
CONLIM: This section sets up various limits 
see comments in Appendix 1 for details 
PROGRAM PROCEDURE: 
discussed in section 3 
PARAMETERS: 
Input for HGEHOG 
CROSS-SECTION DATA: 
PEAR language routine discussed in 
section 4 
S¥YSOUT: Sets up ouput units for various parts 
of the package (see GROUPS in comments) 
OUTFLG: Determines amount of output produced in 
various parts of the package 
(see GROUPS in comments - Appendix 1) 
BODY: Input of SHwave station data and 
observed phases, also precision of data 
P BO: Input of P wave station data and 


observed phases, also precision of data 
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ICNL: Sets the limits for ICNLIM which determines 
how many observation exceed the limits 
imposed by Conlim (see Comments in Blockdata) 

+ GRAV: Input of Gravity data 

LOVE DATA: 

Experimental phase velocities and 
periods or frequencies at which measured 
also precision of data 

RAYLEIGH DATA: 

Experimental phase velocities and 
periods or frequencies at which measured 
also precision of data 

* QDATA: Love and Rayleigh 1/Q data and its precision 

END OF DATA: 

End of this input segment 

SLOP s 


End of program 


The other alternatives are listed in section 2.1 
under the routine INPUT, but a sanien ee such as the above 
covers most of what is commonly useful - of course not 
all these groups have to be present. The details of the 
formats required are to be found in the source program. 


* A feature of KBCAM not implemented at University of 


Alberta 
+ A feature under development 
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